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ABSTRACT 


Several techniques for estimation of the power spectral density (PSD), or 
simply the spectrum, of discretely sampled short data sequences are 
investigated in this thesis. These techniques are directly based on the fast 
Fourier transform (FFT) of the data sequence. However, an extrapolated 
version of the original data sequence is used instead of the original. The aim ‘ 
is to enhance the resolution of the true spectral components of the signal 
without adding any false frequency components. Several data extrapolation 


models are analyzed and their performance at different signal to noise ratios 


(SNRs) and model orders are examined. 
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I. INTRODUCTION 


This thesis investigates several techniques for enhancing the spectrum 
estimation of short data sequences. A short data sequence is defined as one in 
which the spectral resolution required is of the same order, or better, as the 
reciprocal sequence length. The techniques to be investigated will involve 
extrapolating the original data sequence and then performing a basic 
periodogram on this new data sequence. By extending the original data 
sequence we hope to enhance the true spectral components and the spectral 
resolution of the signal without adding any false spectral components. 

The method of spectrum estimation to be used on all extended data 
sequences is the basic periodogram. The periodogram is defined as 


2 
B,(e/®) = “ (1.1) 


M-1 
Dlr 
=0 








where x[n] is the extrapolated data sequence and M is the number of points of 
the extended data sequence. 

This thesis is arranged into 5 chapters and 4 appendices. Chapter II is 
devoted to data extrapolation via linear predictive filtering. In Chapter III the 
data is modeled as the impulse response of an infinite impulse response (IIR) 
filter. Data extrapolation by dominant eigenvectors is discussed in Chapter 
IV. In the last chapter some general conclusions of the work carried out in 
this thesis are presented, and suggestions for future investigations are given. 
Appendix A covers several other methods of AR parameter estimation. 


Appendix B gives some results from a causal Prony's method of 








extrapolation. Appendix C contains some results from some miscellaneous 
test sequences while computer simulation programs are included in 


Appendix D. 





II. LINEAR PREDICTIVE FILTERING 


A. DEFINITION 

In this chapter, all original test data sequences will be extrapolated by 
linear prediction. In linear prediction, a given set of sample values of a 
random sequence x is used to predict the values of the sequence some time 
into the future. We assume that only P previous values of the sequence are 
used in the prediction. Now the problem is to find the prediction coefficients 


—a1, —a9, ..., ~ap to produce an estimate of the current value of the random 


process x[n] from P previous values. The estimate is written as 


£[n] = -a}x[n—1]-a>x[n-2F.. —ayx[n -P] (2.1) 


where * represents the complex conjugation. 


The error in the estimate is given by the difference 
P 
e[n] = x[n]- 2[n] = ¥a,x[n-k] (2.2) 
k=0 


where apo is defined as 1. 
The approach is to choose the prediction filter coefficients to minimize 


the mean-square error 
of = Elletn}?| = E|[+{n] - sin] (2.3) 


To find the optimal prediction filter coefficients we apply the Orthogonality 
Theorem (Ref. 1:p. 308 ] which states that 


E{x[n—iJe*[n]]|=0 fori=1,2,...P (2.4) 





or 











P 
Tyx{i]=—- > ar {i- 4] for i=1,2,...,P (2.5) 
f=1 
where fr, is the autocorrelation function of x. 


The minimum prediction error is found to be 


€min = Txx{0]+ ene (2.6) 
k=] 


These equations are identical to the Yule-Walker equations for autoregressive 
(AR) processes [Ref. 2:p. 157]. Therefore, the optimal linear prediction 
coefficients are just the AR parameters, when the order of the AR process is 
identical to the order of the linear predictor. 


B. AR PARAMETER ESTIMATION 

For good estimates of the AR parameters or prediction coefficients, a 
maximum likelihood estimator (MLE) is usually employed. Several methods 
were tried in this thesis including the autocorrelation method, the covariance 
method, the forward-backward method (modified covariance method), the 
modified forward-backward method, and the Burg method. The two 
methods which proved to be the most robust were the forward-backward 
method and the modified forward-backward method. Explanations and 
results from the other methods are given in Appendix A. 


1. Forward-backward Method (Modified Covariance) 


For an AR process the optimal forward predictor is given by Eq. (2.1). 
The corresponding optimal backward predictor is given by [Ref. 2:p. 226] 


xn] = eee + k}. (2.7) 
k=1 


The forward-backward method estimates the prediction coefficients by 
minimizing the average of the estimated forward and backward prediction 
error powers. The average of the estimated forward and backward prediction 
error powers is given by 








p=5(6' +2) (2.8) 


where the forward and backward prediction error powers are respectively 


defined as 


a ] N-1 


a 


P 2 
x[n]+ ¥ afk]x{n- i 
k=1 





P 2 
NoP x[n]+ Sia [k]x[n +k] , (2.9) 


k=1 








n=0 
and N is the number of data points in the random sequence x. To minimize 
Eq. (2.8) we differentiate 6 with respect to the real and imaginary parts of a[k] 
for k=1,2,...P and set the expression equal to zero. After simplification this 


becomes 


Pp N-1 N-1-P 
Ete 2 2" —k}x*[n-¢]+ Six*[n+k]x[n+ a| 


N-1 oe 
- Laln}e[n— a+ Ee"Inhtne (2.10) 


for 4=1,2,...,P. 
Defining 


; 1 N-1 : : N-1-P Bo 
calitl- aml 3 [n— j]x{n-k]+ Lalns ib neu] (2.11) 


Eq. (2.10) can be written in matrix form as 


CaxlL 1] Corll2] ++ CexltP]] af]]  [Craf1,0] 


Coxl2,1] Cypl2,2] --- Cexl2,P] al2}} __| Cexl2,0] | (2.12) 


CexlP.1] CyalP.2] = CoelP.PIL API] [CeefP-0] 


The first matrix is a P-by-P autocorrelation matrix while the last 
matrix is called the cross correlation vector. The prediction coefficients are 
easily solved for from Eq. (2.12). 

Intuitively, the order P should be as large as possible in order to have 
a large data base for the predictor. The use of too large a value of P can give 
rise to spurious spectral peaks in the AR spectral estimation using the 
modified covariance method [Ref. 3:p. 347]. For the best performance of the 
forward-backward linear predictor (FBLP) method, Lang [Ref. 4:p. 720] suggests 
the value 


P (2.13) 


IR 


N 
3 


where once again N is the original data length. 
Akaike [Ref. 5:pp. 203-217] proposes two other methods of model 
order selection based on the estimated prediction error power. 


a. Test Data 

Two types of test data were used for computer simulation of all 
methods. The first set was obtained from a paper by S. M. Kay and S. L. 
Marple [Ref. 6:p. 1380-1414]. The data is made up of three sinusoids and a 
colored noise process which is obtained by filtering a white Gaussian process. 
The three sinusoids are at fractional frequencies of the sampling frequency of 
0.10, 0.20, and 0.21 and have SNR's of +10, +30, and +30 dB respectively. SNR 
is defined as the ratio of the sinusoid power to the total power in the 
passband noise process. The noise process passband is centered at 0.35 and has 
a bandwidth of 0.30. Figure 1 is a plot of this 64 point test sequence and the 
true PSD. Figure 2 is the periodogram of the original, nonextrapolated data 
sequence using a Hamming data window. A Hamming window is used in all 
of the periodogram estimations in this thesis. 
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Figure 1. Kay and Marple Data Set with True Power Spectral Density 
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Figure 2. Periodogram of Kay and Marple Data (Hamming Window) 





The second type of test data is made up of two sinusoids 
corrupted by white noise . This is also a 64 point test sequence. This type of 
test data will be composed of different frequencies, phases, and SNRs to see 
how each method of spectrum estimation resolves the sinusoidal frequencies. 
A sampling frequency of 1000 Hz is used to generate the signals, therefore the 
maximum frequency of any single sinusoid is 500 Hz. Using ordinary FFT 
routines on non-extended 64 point data sequences, the spectral resolution can 
be as good as 0.8(1/64msec)=12.5 Hz. Using extended data sequences, it will be 
shown that the resolution can be improved to 2.0 Hz for an appropriate SNR. 

For comparison purposes, the periodograms of the extrapolated 
data sequences are plotted with the autoregressive (AR) spectrum estimations 
of the same data sequences. The AR spectrum is defined as 


Sar(e/#) = —1 (2.13) 
where 
; P 3 
Aei@) = a{njeio (2.14) 
n=0 


and of course the AR coefficients are identical to the prediction filter coefficients. 
The dashed lines on each plot represent the positions of the true 
frequencies of the test data. 
b. FBLP Results 
Figure 3 is the extrapolated version of the Kay and Marple (KM) 
data. The data was extrapolated by using linear prediction as in Eq. (2.1). The 
prediction coefficients were obtained from the FBLP method utilizing Eq. 
(2.12). The data length was extrapolated to 320 points using a model order of 
16. 











Care must be taken when extrapolating the data sequence. If the 


correlation matrix in the calculation of the prediction coefficients is not 
positive definite, the resulting poles in the filter of the AR model are not 
guaranteed to lie within the unit circle. This could cause the extrapolation to 
grow without bound and will have devastating effects on the spectrum 
estimation. The forward-backward method does not guarantee a stable all 
pole filter, although in practice the extrapolation is usually bounded. We 
should check to ensure all poles of the AR process lie within the unit circle 
before we extrapolate the data. 

Figure 4 shows the periodogram of the extrapolated KM data and 
the respective AR spectra] estimation. Both methods correctly calculate the 
three sinusoidal frequencies present. The colored noise process is also given a 
fair representaion. When the data is extrapolated to 576 points, an even better 
PSD estimate is obtained from the FBLP extrapolation. Figure 5 shows the 
results. 

For the two sinusoids corrupted with white noise, good results 
are obtained. Figure 6 shows the spectrums of two sinusoids with a 10 Hz 
separation and a SNR of 12 dB. The original data sequence was extrapolated 
to 320 points. Both sine waves have a 0 degree initial phase angle. The 
frequencies are correctly calculated using the periodogram of the FBLP 
extrapolated sequence while the AR spectral estimation is unable to 
distinguish the two closely spaced sinusoids. 

When the initial phase angle between the two sine waves is 
changed to 45 degrees, the periodogram of the FBLP extrapolated sequence, 
Figure 7, distinguishes the two frequencies, even though one of the estimated 
spectral components is slightly off. The AR spectrum estimation is unable to 
resolve them. Once again the data was extrapolated to 320 points. 
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Figure 3. Extrapolated KM Data, Model Order = 16 
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Figure 4. Periodogram of Extrapolated KM Data (320 points) and Respective 
AR Spectrum Estimation 
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Figure 5. Periodogram of Extrapolated KM Data (576 points) and Respective 
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MODIFIED COV EXTENSION, ORDER 16, SNR=12, PIT=0 
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Figure 6. Periodogram of Extrapolated Sequence and Respective AR 
Spectrum Estimation of two Sine Waves (0° Initial Phase Angle) 
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Figure 7. Periodogram of Extrapolated Sequence and Respective AR 
Spectrum Estimation of two Sine Waves (45° Initial Phase Angle) 
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When the initial phase angle between the two sine waves is 

changed to 90 degrees, in essence a signal consisting of a sine wave and a 

cosine wave with frequencies that are 10 Hz apart, neither method is able to 
distinguish the two frequencies. 

2. Modified Forward-Backward Method 

The modified forward-backward method is a two-step procedure for 

modifying the FBLP method [Ref. 3:p. 347]. First, we compute the P 

eigenvalues of the correlation matrix given in Eq. (2.12). To do this, we use 


the singular value decomposition (SVD) which is applied to the 2(N-P)-by-P 


data matrix A given by 
x{P] x[P+1] --- x{[N-1] x*[2] x*[3] + x*[N-P+]] 
AH +P -1} +P] se x[N -2] x*[3] x*[4) “+ x*[N ~P +2] 
x(t] x(2] 5 x{N-P] x*[P ay] x*[P+2] oe x*{N] 


(2.15) 
where N is the length of the original data sequence. Note that the correlation 
matrix in Eq. (2.12) is equal to AHA. 

Once the P eigenvalues are computed, we divide the data space into 
two parts: (i) the signal subspace spanned by the eigenvectors associated with 
the K largest eigenvalues, and (ii) the noise subspace spanned by the 
remaining (P-K) eigenvectors. The notion used is that the K principal 
eigenvectors of the correlation matrix are due to the signal components and 
are less sensitive to perturbations caused by noise than the remaining 
eigenvectors. In the case of noiseless data consisting of K sinusoids, the 
correlation matrix has K nonzero eigenvalues and (P-K) zero eigenvalues. 
Thus in a practical case of noisy data, by keeping the K largest eigenvalues of 
the correlation matrix and ignoring the rest, we are in effect attempting to 


increase the SNR by approaching the noiseless case. 
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Let v1, 02, ..., Vx denote the K largest eigenvalues and e1, €2, ..., ex 


denote their associated eigenvectors, respectively. We use these eigenvalues 
and eigenvectors to compute the following P-by-1 estimate of the prediction 
coefficients: 
K éi 
a= >i“Lef Mp) (2.16) 
i=1 Vi 
where 
b =[x[P+1],x[P +2],...,x[N],x *[1],x*[2],....x*[N —P]]. (2.17) 


For the modified FBLP (MFBLP) method, we have removed the 
undesirable effects of the noise subspace eigenvectors, therefore we may 
increase the predictor order P. Tufts and Kumaresan [Ref. 7:pp. 975-989] have 
experimentally determined the value of the predictor order to be 


P=3N/4 (2.17) 


where N is the original, nonextrapolated data length. 

Since the computer simulation uses 64 data points, a model order of 
48 is selected. 

a. Modified FBLP Results 

Many test runs are made on the signal consisting of 2 sine waves 

corrupted by white noise. Sinusoidal frequencies, phases, and SNRs are 
changed for comparison purposes. Table 1 specifies the parameters for each 
run and the corresponding figure number. In each test case, the data sequence 
is extrapolated to a data length of 320 points using linear prediction with the 
prediction coefficients obtained from the MFBLP method. Figures 8 through 
24 show the results. 
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TABLE 1. PARAMETERS FOR TEST SEQUENCES 
CASE #{ SINUSOIDAL SNR (DB) RELATIVE 
FREQUENCIES 
(HZ) 


INITIAL PHASE 
| 1 | 100and107_—| 2 
| 2 | 100and 107 _| 
| 3 | 100and 107_| 
| 4 | 100andi07_—| 
| 5 | 100 and 107_| 
| 6 | 100and 107__| 
| 7 | 100and 107_| 
| 8 | 100and107_ | 
a 
| 10 
ae 
| 12 
paige) 
ae 
p15 





FIGURE # 
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sequence performs as well or better than the AR spectrum estimation. Figure 
8 shows that even though both methods correctly estimate the frequencies, 
the periodogram of the MFBLP extrapolated sequence sometimes has better 
defined spectral peaks. In case 10, depicted in Figure 17, the periodogram of 
the MFBLP extrapolated sequence correctly distinguishes the two frequencies 
while the AR spectrum estimation is unable to do so. Even in test case 17 
where the two sine waves are only 2 Hz apart, the periodogram of the MFBLP 
extrapolated sequence has well defined peaks, as shown in Figure 25, while 


the AR spectrum estimation just barely distinguishes the two frequencies. 
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Figure 8. Case 1 Spectrum Estimates 
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Figure 9. Case 2 Spectrum Estimates 
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MODIFIED FBLP EXTENSION, ORDER 48, SNR=6, PH=0 
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Figure 10. Case 3 Spectrum Estimates 
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Figure 11. Case 4 Spectrum Estimates 
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MODIFIED FBLP EXTENSION, ORDER 48, SNR=9, PI1=45 
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Figure 12. Case 5 Spectrum Estimates 
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MODIFIED FBLP EXTENSION, ORDER 48, SNR=6, PH=45 
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Figure 13. Case 6 Spectrum Estimates 
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MODIFIED FBLP EXTENSION, ORDER 48, SNR=3, PH=45 
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Figure 14. Case 7 Spectrum Estimates 
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Figure 15. Case 8 Spectrum Estimates 
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MODIFIED FBLP EX TENSION, ORDER 48, SNR=9, PIT=90 
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Figure 16. Case 9 Spectrum Estimates 
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Figure 17. Case 10 Spectrum Estimates 
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Figure 18. Case 11 Spectrum Estimates 
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FBLP EXTENSION, ORDER 48, SNR=12, PH=0 
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Figure 19. Case 12 Spectrum Estimates 
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MODIFIED FBLP EXTENSION, ORDER 48, SNR=12, PH=45 
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Figure 20. Case 13 Spectrum Estimates 
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Figure 21. Case 14 Spectrum Estimates 
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MODIFIED FBLP EXTENSION, ORDER 48, SNR=12, PH=90 


a 
acd 
a) 
na 
a 
& 
< 
—t 
iy) 
a 


0.15 0.2 0.25 0.3 0.35 0.4 


FRACTION OF SAMPLING FREQUENCY 


AR SPECTRUM ESTIMATION 


RELATIVE PSD (dB) 


0.15 0.2 0.25 0.3 0.35 U.4 0.45 
FRACTION OF SAMPLING FREQUENCY 





Figure 22, Case 15 Spectrum Estimates 
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Figure 23. Case 16 Spectrum Estimates 
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Figure 24. Case 17 Spectrum Estimates 
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Figure 25. Periodogram of Extrapolated Data Sequence for Case 17 





When the MFBLP extrapolation technique is tested on the KM 


sequence, the periodogram of the extrapolated sequence is considerably better 
than the AR spectral based estimation. In general, the modified FBLP method 
utilizes prior knowledge of the number of sinusoids present in the signal. 
When the K largest eigenvalues are chosen, we are assuming that there are K 
sinusoids present in the signal. If one of these K largest eigenvalues is 
associated with the noise, inaccurate estimates will be computed. This can be 
seen in Figure 26. Because the KM data is made up of 3 sinusoids, the 3 
largest eigenvalues and corresponding eigenvectors are chosen. However in 
this case, one of these eigenvalues is associated with the noise so the signal at 
the 0.1 fractional frequency is not detected in the AR spectral estimation. 
However when the same prediction coefficients are used in the data 
extrapolation, the resulting periodogram is a more accurate spectral 
estimation than the AR spectral estimation. Obviously, the extrapolation 
approach is more forgiving if the exact number of sinusoids is unknown. In 
Figure 27 the four largest eigenvalues and corresponding eigenvectors were 
chosen, one more than the actual number of sinusoids present. Once again 
the AR spectral estimation is incorrect, picking up more noise, while the 
periodogram of the modified FBLP extrapolation gives a reasonable spectral 
estimate. The KM sequence is extrapolated to a length of 320 points before the 


periodograms are computed. 
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MODIFIED FBLP EXTENSION, ORDER 48, SNR=30, PH=0 
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Figure 26. KM Spectrum Estimates Using the Three Dominant Eigenvectors 
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Figure 27. KM Spectrum Estimates Using the Four Dominant Eigenvectors 
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Ill. PRONY'S METHOD 


A. DEFINITION 
The second method of sequence extrapolation involves modeling the 
data as the impulse response of a linear time invariant filter. This approach 


is shown in Figure 28. Note that B(z) has q zeros and A(z) has p poles. 





Figure 28. Impulse Response Estimation 


We want to minimize the error e[n] in some sense. If we would 
minimize Sleny’, then this leads to non-linear equations and is very 
difficult to solve. Therefore we will use an indirect method to solve this 


problem. The equation for the filter is 
X[n]+ a,X[n-1]+...+a,x[n — p] = by 5(n) + b)5(n—1)+...+b,5(n—q). (3-1) 


The coefficients aj and bj are unknown in Eq. (3.1) and we want the output 


estimate x{n] as close to x{n] as possible. 


Assume that we have N data points and that N 2 p+q+1. The transfer 
function H(z) of the filter is 


H(z) = Biz) _ +byz'+...+b,274 BS 
A(z) 1+a27}+.. yz P 


We can write Eq. (3.1) in matrix form as 


x{0] 0 O- 0 Po 
fo} 0 -- 0 : ; 
Ag a-t] fa-p]_ [|=], (3.3) 
x[q+1] xg) oe x[qt+1-p} } : 4 
: Qy : 
x[N-1] x[N-2] -- + x[N-1-p] 


or more compactly 
XB b 
eaPele os 


For the case where N=p+q+1, we can write the lower set of equations as 


where 0° = [0, 0,..., 0]. 


X,a=0 (3.5) 


and solve directly for a. Once we know a we can solve the upper equations 
directly for b where 
Xpa = b. (3.6) 


This method is known as Pade's approximation [Ref 1:p. 495]. 
For the more general case where N > p+q+1 we still have Eqs. (3.5) and 
(3.6), however now Eq. (3.5) is an overdetermined set. If we solve Eq. (3.5) in a 


least squares sense, then we can still find b from Eq. (3.6). Solution of this 
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least squares problem [Ref. 1:pp. 474-495] leads to the least squares normal 


equations 


(x4 X,)a= | (3.8) 


where S = Y, len]. 


This is known as Prony's method for time domain IIR filter design. This 
method was used to model the test data as the impulse response of an IIR 
filter. This impulse response is then extrapolated over a longer duration and 
then its periodogram is computed. 


B. NONCAUSAL EXTRAPOLATION 

Figure 29 shows the original Kay and Marple data sequence and the 64 
point Prony's impulse response estimate of this data sequence. The order of 
the numerator and denominator in the transfer function of Eq. (3.2) was 
chosen to be 16. Several model orders were tried before settling on the order 
of 16. If too low an order is chosen, there is not enough resolution in the 
periodogram of the extrapolated sequence. If too high an order is chosen, 
spurious spectral peaks arise. Examples of varying model orders are given in 
Appendix B. As can be seen in Figure 29, Prony's impulse response estimate 
of the data sequence for an order of 16 is a very good approximation. 

By choosing the appropriate impulse duration, the impulse response can 
be extrapolated as far as we like. The impulse response can also be 
extrapolated in both forward and backward directions. We can extend it 
backwards by reversing the original data sequence and performing the same 
computations as one normally would for the forward case. Once this new 
Prony's estimate has been computed, we reverse the extension and attach it to 
the front of the original data sequence. This gives a noncausal extrapolation 


of the original data sequence. 
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Figure 29. Original KM Data and Prony’s Impulse Response Estimate of KM 
Data 
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This was the method performed on all test sequences. A strictly causal or 
forward extrapolation was tried but the results were not as good as the 
noncausal extrapolation. Some examples of the causa] extrapolation are 
shown in Appendix B. The length of the extrapolation is provided on each 
figure. 


C. RESULTS 

For the KM data, the noncausal extrapolation method gave excellent 
results. Figure 30 shows the extrapolated time series and the resulting 
periodogram. The three sinusoids present are correctly identified and the 
colored noise passband process is properly represented. 

Figure 31 shows the extrapolated sequence and corresponding 
periodogram of a test sequence consisting of two sinusoids embedded in 
white noise. The two sinusoids have frequencies which are 10 Hz apart 
(fs; = 1000 Hz) and both have zero initial phase. The SNR is 15 dB. The 
periodogram of the extrapolated sequence correctly identifies the two 
frequencies present. When the initial phase of one of the sine waves is 
changed to 45 degrees, the periodogram of the extrapolated sequence once 
again identifies two frequencies, however one of the frequency estimates is 
slightly off. Figure 32 shows the results. Figure 33 shows the method does 
not distinguish the two frequencies when one of the initial phase angles is 
changed to 90 degrees. 

Figures 34 through 36 show the extrapolated time series and resulting 
periodograms for two sine waves 5 Hz apart in frequency with different initial 
phase angles. The SNR is 30 dB in all three cases. The periodograms in each 
of the three cases correctly distinguished the two frequencies present when 
the initial relative phases were 0 degrees (Figure 34), 45 degrees (Figure 35), 
and 90 degrees (Figure 36). 


PRONYS EXTRAPOLATED ESTIMATE - ORDER 16 NON-CAUSAL 


amplitude 


160 180 200 


sample 


PRONYS METIIOD PERIODOGRAM - ORDER 16 NON-CAUSAL 


4 
Han 
fe 
hfe 
te 
eye 
’ 


. 
fe 
ode 
ve 


oe 
te 
oe 
16 


oo 
ot 
et 


magnitude in dB 


0.1 O.15 0.2 0.25 0.3 0.35 


fraction of sampling frequency 





Figure 30. Extrapolated KM Data and Respective Periodogram 
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Figure 31. Extrapolated Sequence of Two Sine Waves, 10 Hz Apart, 0° Initial 
Phase Angle, SNR = 15 dB, and Resulting Periodogram 
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Figure 32. Extrapolated Sequence of two Sine Waves, 10 Hz Apart, 45° Initial 
Phase Angle, SNR = 15 dB, and Resulting Periodogram 
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Figure 33. Extrapolated Sequence of two Sine Waves, 10 Hz Apart, 90° Initial 
Phase Angle, SNR = 15 dB, and Resulting Periodogram 
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Figure 34. Extrapolated Sequence of two Sine Waves, 5 Hz Apart, 0° Initial 
Phase Angle, SNR = 30 dB, and Resulting Periodogram 
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Figure 35. Extrapolated Sequence of two Sine Waves, 5 Hz Apart, 45° Initial 
Phase Angle, SNR = 30 dB, and Resulting Periodogram 
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Figure 36. Extrapolated Sequence of two Sine Waves, 5 Hz Apart, 90° Initial 
Phase Angle, SNR = 30 dB, and Resulting Periodogram 
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The periodogram resulting from Prony's method of data extrapolation is 
unable to distinguish the two sinusoids which are 5 Hz apart in a SNR less 
than approximately 27 dB. Prony’s method of data extrapolation resulted in 
periodograms which were not as accurate as the MFBLP method in lower 
SNRs. In high SNRs (20 dB or higher) good results are obtained but the 
periodograms resulting from the data extrapolation using the MFBLP 


prediction coefficients still provide better estimates. 
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IV. EIGENVECTOR EXTRAPOLATION 


A. DEFINITION 

In this chapter we develop a series expansion for a random vector in 
terms of an arbitrary set of orthonormal vectors [Ref. 8:pp. 261-262]. Once this 
set of orthonormal vectors (eigenvectors) are computed, we will use the 
dominant ones to extrapolate the data sequence and, hopefully, improve the 
spectral resolution. 

First consider a set of N orthonormal vectors qu, q2, ..., qn that satisfy the 
following two requirements: 


Lok=% 
T, Jb 


Let x denote an N-by-1 random vector with zero mean and correlation matrix 
R,. The requirement is to express the random vector x as a linear 


combination of the set of orthonormal vectors qui, q2, ..., N as follows: 
N 
x= Y aai- (4.2) 
i=] 
Equations (4.1) and (4.2) define the discrete time version of the Karhunen- 
Loéve expansion [Ref. 8:p. 261]. The coefficients of the expansion, aj, are 
defined by the inner product between the vectors x and qj, given by 
a = qi x 
=x'q; fori=1,2,...,N. (4.3) 
The coefficients a; are obviously random variables and we would like to 


choose the qj so that these coefficients are mutually uncorrelated. 
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Because the random vector x has zero mean, we observe that all the 


coefficients a; also have a mean of zero, that is 
E[a;]=0 fori=1,2,...,N. (4.4) 


Therefore, for the coefficients of expansion to be uncorrelated, we require that 


they satisfy the condition 


cov[a,.a;] = E|(o, — Ela, ](a; - E{a;})| 
= E[a,0%;] 


_ | Uj, k=i 
{0 k#i 9) 


where the vj; are constants to be specified. From Eq. 4.3 we can rewrite Eq. 4.5 


as 
cov| a%,0,;| = E[a,e;| 

= Elapxx"q;| 

= aE E[xx" Jai (4.6) 
Since by definition the correlation matrix of the random vector x equals 

R,= E|xx™ |, (4.7) 

we can simplify Eq. 4.6 as 

cov|a,a;] = qf Rqi- (4.8) 
Therefore from Eq. 4.5 and Eq. 4.8 we have 


v;, k =i, 
afRai-{") ‘ee (4.9) 


If we multiply both sides of Eq. 4.1 by v;, we have 





vj, k=i, 
ViGkGi -{ Ses (4.10) 
Therefore for Eq. 4.9 to hold for all choices of k and a particular i, it is 
necessary that the matrix product R,qi equal v;qi, as shown by 


R,.qj = 0;4; (4.11) 


where the scalar v; is the eigenvalue of the correlation matrix R, and the 
vector qi is the eigenvector associated with the eigenvalue 1}. 

Therefore we can express the random vector x as a linear combination of 
the eigenvectors of the correlation matrix of x. Once we have computed the 
eigenvalues and eigenvectors, we will extract the dominant ones and use 
them to extend the original data sequence by appending this new sequence to 
the original sequence. All test sequences are extrapolated to a length of 128 


points. 


B. RESULTS 

The first set of test data consists of two sinusoids with a 10 Hz (f, = 1000 
Hz) difference in frequency and 0 degree relative phase angle. By examining 
the eigenvalues of the correlation matrix, we observe that there are two 
dominant eigenvectors. This corresponds to the fact that we have two 
sinusoids present. Because there are two dominant eigenvectors present, we 
use the linear combination of the two most dominant eigenvectors for the 
extrapolation. Therefore in Eq. 4.2, N is set to 2 and the resulting sequence is 
appended to the original data. The autocorrelation matrix in Eq. 4.2 is 
computed by using the biased autocorrelation method. Figure 37 shows the 
plot of the extrapolated data sequence and the resulting periodogram. Two 
distinct peaks are present but they are slightly off. The SNR is 15 dB. 
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Figure 37. Eigenvector Extrapolation of 2 Sine Waves, 10 Hz Apart, 0° Initial 
Phase Angle, SNR = 15 dB, 2 Eigenvectors Used; and Resulting Periodogram 
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Figure 38 shows the results when the initial relative phase angle between 
the two sine waves is changed to 45 degrees. Once again two distinct spectral 
peaks are present but they are slightly off from the true values. Figure 39 
shows the results when the relative phase angle is changed to 90 degrees. As 
in the two prior cases, the spectral peaks are slightly off from their true 
values. 

Even though the spectral peaks are not absolutely accurate, they remain 
fairly consistent in even lower SNRs. Figure 40 shows the results of the two 
sine waves (initially in phase, 10 Hz apart) in a SNR of -3 dB. The spectral 
peaks are still well defined. When the SNR was lowered to -6 dB, the 
periodogram was unable to distinguish the two peaks. Figure 41 shows these 
results. 

Figure 42 shows the results of the eigenvector extrapolation on the KM 
sequence. N was chosen to be 3 as we have three sinusoids present in this 
case. Three main spectral peaks are present but one of the estimates is slightly 


off. 
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Figure 38. Eigenvector Extrapolation of Two Sine Waves, 10Hz Apart, 45° 
Initial Phase Angle, SNR = 15 dB, 2 Eigenvectors Used; and Resulting 
Periodogram 
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Figure 39. Eigenvector Extrapolation of 2 Sine Waves, 10 Hz Apart, 90° Initial 
Phase Angle, SNR = 15 dB, 2 Eigenvectors used; and Resulting Periodogram 
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Figure 40. Eigenvector Extrapolation of 2 Sine Waves, 10 Hz Apart, 0° Initial 
Phase Angle, SNR = -3 dB, and Resulting Periodogram 
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Figure 41. Eigenvector Extrapolation of 2 Sine Waves, 10 Hz Apart, 0° Initial 
Phase Angle, SNR = -6 dB, and Resulting Periodogram 
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Figure 42. Eigenvector Extrapolation of KM Data, 3 Eigenvectors Used, and 
Resulting Periodogram 
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V. CONCLUSIONS AND RECOMMENDATIONS 


In this thesis, data extrapolation techniques to enhance the spectral 
resolution of the periodogram are presented. The three methods of 
extrapolation are: linear predictive filtering, modeling the data as the 
impulse response of a filter (Prony's method), and eigenvector extrapolation 
using the Karhunen-Loéve expansion. 

The most promising of the three methods is the linear predictive filtering 
extrapolation using the prediction coefficients obtained from the modified 
forward-backward method. This technique resulted in spectral estimates 
which were equivalent to or better than those obtained via the autoregressive 
method. Frequencies as closely spaced as 2 Hz have been resolved, with an 
appropriate SNR, where normally resolution of no better than 12.5 Hz would 
be obtained using a periodogram of the nonextrapolated sequence. This 
method performed very well in low SNRs. 

In a high SNR environment, the periodograms resulting from Prony’s 
method of extrapolation gave good results. Performance was not as 
promising in lower SNRs. Model order selection in Prony’s method proved 
to be very important in determining the resolution of the resulting 
periodogram. 

None of the methods covered required very long computer run times. A 
typical run lasted approximately 1 to 2 minutes on a 286 based personal com- 
puter, depending on the extrapolation length. As discussed earlier, care must 
be taken when extrapolating the data sequence to ensure that it does not grow 


without bound. This has devastating effects on the spectrum estimation. 
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In this thesis we have shown that data extrapolation is a viable means for 


increasing the resolution of the periodogram. Further research should 
concentrate on other data extrapolation techniques for better and more robust 


spectral resolution. 








APPENDIX A. OTHER METHODS OF AR PARAMETER ESTIMATION 


A. AUTOCORRELATION METHOD 
In the autocorrelation method or the Yule-Walker approach, the AR 
parameters are estimated by minimizing an estimate of the prediction error 


power 


, 1 
pea 


n=—eo 


Pp 2 
x[n]+ dalk}e{n —k] . (A.1) 








The samples of x[n] which are not observed, those not in the range 0Osn<N-1, 
are set equal to zero. To minimize the prediction error power, we 
differentiate Eq. A.1 with respect to the real and imaginary parts of the alk]s. 
This yields 

| aes P 

a Y [sine Eaeain-n et -—f=0, 4£=1,2,...P. (A.2) 


n=—0o =1 
In matrix form this set of equations becomes 


Tx[0] Txx{-1] ad Tzx[-(P ~1)] a{1] rxx{1] 


ralll lO) ~~ real-(P~2i)] 421] |ral2l] gy 
relP—1] telP-2)-~ reall fala]] [rel PI 
where 
1 N-1-k 
rel] =1N 2e"[h{n+#] for k = 0,1,...,P (A.4) 
rel fork =-(p—1),-(p-2)ponr-1 
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which is recognized as the biased autocorrelation function (ACF) estimator. 
The AR parameters or prediction coefficients are easily solved for from Eq. 
A.3. 

The matrix in Eq. A.3 is hermitian (rxx{-k] = r,[k]) and Toeplitz, and 
furthermore can be shown to be positive definite. As such, the resulting 
estimated poles are guaranteed to be within the unit circle. This property 
ensures that the prediction coefficients calculated from Eq. A.3 will give stable 
extrapolations. This is very important because if the extrapolation grows 
without bound, it will have devastating effects on the spectrum estimation. 

The autocorrelation method extrapolation was found to produce poorer 
resolution of the spectral estimates than the other extrapolation techniques 
used. The length of the extrapolated sequences in this section is 320 points. 


Figure 43 shows the resulting periodogram of the autocorrelation method 


extrapolation and the respective AR spectrum estimation of the KM data. 


The extrapolation gave a periodogram which correctly identifies the 3 
sinewaves present but the spectral peaks are not as well defined as in the 
periodogram obtained from the MFBLP extrapolation technique. The AR 
spectrum estimation resulting from the autocorrelation method is unable to 
distinguish the two closely spaced sinusoids at 0.2 and 0.21. 

When the autocorrelation method extrapolation was applied to a data 
sequence consisting of two sine waves embedded in white noise (SNR=15 dB), 
the resulting periodogram was unable to distinguish the two closely spaced 
sinusoids (5 Hz apart, fs = 1000 Hz). The AR spectrum estimation also failed. 


Figure 44 shows the results. 
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Figure 43. Periodogram of Extrapolated Sequences and AR Spectral 


Estimation of KM Data Using the Autocorrelation Method 
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Figure 44. Periodogram Resulting from Extrapolation of 2 Sine Waves, 5 Hz 
Apart, SNR = 15 dB, 0° Initial Phase Angle and Respective AR Spectrum 





Estimation 


68 








B. COVARIANCE METHOD 


In the covariance method, the AR parameters are estimated by 


minimizing the estimate of the prediction error power 





5 Slates Sateen (a5) 
=-—_—_—— xin\+ a xin— é . 
: N-P Jp k=1 





Note that the only difference between the covariance method and the 
autocorrelation method is the range of summation in the prediction error 


power estimate. The minimization of Eq. A.5 yields 


Cxx{11] yy [1,2] = cxe{1,P] | aft] Cxx[1,0] 


ceL2A} coal22] ~~ Cel2-PH] Zl] Jem] ig g 
cell] CeelP.2] -~ CgelPPIL PI] [eal PO 
where 
1 N-1 
cali Alig Sette bln (47) 


The correlation matrix in Eq. A.6 is positive semidefinite and not 
Toeplitz. As such, the estimated poles using the covariance method are not 
guaranteed to lie within the unit circle. Therefore the prediction coefficients 
are not guaranteed to give bounded extrapolations although in practice the 
extrapolations are generally bounded. 

Results from the covariance method extrapolation were better than the 
autocorrelation method extrapolation. Once again the data is extrapolated to 
a length of 320 points. Figure 45 shows the periodogram of the covariance 
method extrapolation and the respective AR spectrum estimation of the KM 
data. The periodogram from the covariance method extrapolated data se- 
quence has well defined spectral peaks and correctly identifies the 3 sinusoids 
presen: The AR spectrum estimation also identifies the 3 sinusoids. 
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COVARIANCE METHOD - EXTRAPOLATION, ORDER 16 
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Figure 45. Periodogram of Covariance Method Extrapolated KM Data and 
Respective AR Spectrum Estimate 
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When the covariance method extrapolation was applied to a data 
sequence consisting of 2 sine waves embedded in white noise (SNR=15 dB), 
the resulting periodogram was unable to distinguish the two closely spaced 
sinusoids (5 Hz apart). The AR spectrum estimation also failed. Figure 46 


shows the results. 


C BURG'S METHOD 

In contrast to the autocorrelation and covariance methods which estimate 
the AR parameters directly, the Burg method estimates the reflection 
coefficients and then uses the Levinson recursion to obtain the AR parameter 
estimates. 

Burg's optimization criterion is to minimize both the forward prediction 


error ep(n) and the backward prediction error ep(n): 
N-1 
1 +7..72 —7.72 
= —— ) |ep[n]*+ep[n (A.8 
nap Lleol? +r" ) 


where 


ep(n] = x[n] + a[1]x[n — 1] + a[2]x[n- 2]+...+a[P]x[n — P] 
ep[n] = x[n— P]+ af1]x[n- P+ 1]+a[2]x[n— P+ 2]+...+a[P]x]n]. 


The computional steps are summarized below [Ref. 2:pp. 228-231] 
1. Initialize by setting 


eg[n]=eo[n]=x[n] forO<n<N-1 


and the 0th order mean square error as 


N-1 
Eg =~ Dalry’ (A.9) 


n=0 
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COVARIANCE METHOD - EXTRAPOLATION, ORDER 16 
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Figure 46. Periodogram of Extrapolated Sequence and Respective AR 
Spectrum Estimate Using Covariance Method 
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2. Compute the reflection coefficient: 
N-1 7 
2 dep-ilnJep-iln— 1] 
Yp= ipa (A.10) 
Deb-aln} + ep_y[n— 1] 
n=P 


3. Compute ap[k] using the Levinson recursion given by: 


1 0 


1 
apf] ap_,[1] ap_;[P—1] 
ap[2] ~ ori -Yp cae = 2) (A.11) 
ie P] ap_,[P—1] ap_,(1] 


0 J [ 1 


4. Compute ep[n] and ep[n] for P< n <N-1 


ep[n] = ep_s{n]- ypep_[n] 
ep[n] = ep_y[n— 1]- ypep_,[n]. (A.12) 


5. Update the mean-square error as follows 


Ep =(1- yp) Ep-1. 
6. Continue the iteration until P equals the model order. 
Figure 47 shows the results using Burg's method. The data is extrapolated 
to 320 points. The periodogram of the extrapolated KM data resolves the 3 
sinusoidal components, however one of the spectral components is slightly 


off. In addition it seems to be prone to line splitting. The respective AR 


spectral estimation has similar results. The model order is 16. 
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Figure 47. Periodogram of Extrapolated Sequence and Respective AR 


Spectrum Estimate of KM Data Using Burg’s Method 








Figure 48 shows the PSDs of a sequence consisting of 2 sine waves 
embedded in white noise (SNR=15 dB) which are 5 Hz apart (fs =1000 Hz). 
Neither the periodogram of the extrapolated sequence nor the AR spectrum 
estimation could distinguish the two closely spaced sinusoids. The model 


order once again is 16 and the data is extrapolated to 320 points. 
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Figure 48. Periodogram of Extrapolated Sequence and Respective AR 
Spectrum Estimate of 2 Sine Waves using Burg’s Method 








APPENDIX B. SIMULATION RESULTS FROM PRONY'S METHOD 


In Chapter III, the model order for the numerator and denominator in Eq. 
3.2 was chosen to be 16. If too low an order is chosen, there is not enough 
resolution in the periodogram of the resulting extrapolated sequence. Figure 
49 shows the results of a noncausal extrapolation and the resulting 
periodogram for the KM data when the model order is 8. The spectral 
component at 0.1 fraction of the sampling frequency is not resolved and the 
two closely spaced components at 0.2 and 0.21 are incorrect. The data 
extrapolation length is noted on each figure. 

If the order is too high, spurious spectral peaks arise. This can be seen in 
Figure 50 which is the noncausal extrapolation and resulting periodogram of 
two sine waves embedded in white noise (SNR=15 dB) which are 10 Hz apart 
(fs = 1000 Hz). The two frequency components are properly represented but 
an incorrect spectral component at approximately 0.25 arises. 

When a causal or strictly forward extrapolation was tried, the results were 
not as promising as the noncausal extrapolation. Figure 51 shows the results 
of a causal extrapoiation and the resulting periodogram of the KM data. The 
model order is 16. The periodogram of the causal extrapolation was unable to 
resolve the two closely spaced sinusoids at 0.2 and 0.21. Figure 52 shows the 
results of a causal extrapolation and the resulting periodogram of two sine 
waves embedded in white noise (SNR=15 dB) and 10 Hz apart (fs = 1000 Hz). 
The periodogram from the causal extrapolation was unable to resolve the two 


frequencies. 


PRONYS ESTIMATE OF DATA SEQUENCE - ORDER 8 NON-CAUSAL 


3 
4 
= 
= 
E 
o 


PRONYS METHOD - ORDER 8 NON-CAUSAL 


magnitude in dB 


0.2 0.25 0.3 0.35 0.4 0.45 





fraction of sampling frequency 


Figure 49. Noncausal Extrapolation of KM Data and Resulting Periodogram 
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Figure 50. Noncausal Extrapolation and Resulting Periodogram of 2 Sine 
Waves 
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PRONYS ESTIMATE OF DATA SEQUENCE - ORDER 16 CAUSAL 
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Figure 51. Causal Extrapolation and Resulting Periodogram of KM Data 
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Figure 52. Causal Extrapolation and Resulting Periodogram of two Sine 
Waves 
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APPENDIX C. RESULTS FROM MISCELLANEOUS SEQUENCES 


The four main types of extrapolation techniques in this thesis (modified 
covariance, modified forward-backward, Prony, and eigenvector method) 
were tried on a noise only sequence. Figures 53 through 56 show the 
respective results of the periodogram of the extrapolated sequences and the 
AR spectrum estimations. The noise sequences are extrapolated to 320 points 
prior to the calculation of the periodogram. Figure 57 shows the periodogram 
of the nonextrapolated noise sequence. 

The MFBLP method and Prony’s method were also tried on a test 
sequence consisting of a single exponentially decaying sine wave with an 
SNR of 15 dB. The sine wave decayed to 10% of its maximum value within 
the original 64 points. Figure 58 shows the resulting periodogram from the 
MFBLP method extrapolated sequence and the respective AR spectrum 
estimation. Both methods correctly identified the frequency of the sine wave. 
Figure 59 shows the Prony’s method extrapolated sequence and the resulting 
periodogram. Once again the periodogram correctly identified the frequency 


present. 
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Figure 53. Periodogram of FBLP Extrapolation of Noise and Respective AR 
Spectrum Estimate 


MODIFIED FBLP EXTENSION, ORDER 48, NOISE SEQUENCE ONLY 


RELATIVE PSD (dB) 
R 


0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 
FRACTION OF SAMPLING FREQUENCY 


AR SPECTRUM ESTIMATION 


~r 


RELATIVE PSD (dB) 





- a 
9 0.05 0.1 0.15 0. 0.25 0.3 0.35 0.4 0.45 0.5 
FRACTION OF SAMPLING FREQUENCY 


Figure 54. Periodogram of MFBLP Extrapolation of Noise and Respective AR 
Spectrum Estimate 
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Figure 55. Prony’s Extrapolation of Noise and Resulting Periodogram 
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Figure 56. Eigenvector Extrapolation of Noise and Resulting Periodogram 
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Figure 57, Periodogram of the Nonextrapolated Noise Sequence 
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Figure 58. Periodogram of MFBLP Extrapolated Exponentially Decaying Sine 
Wave and Respective AR Spectrum Estimation 
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Figure 59. Prony’s Extrapolation of Exponentially Decaying Sine Wave and 
Respective AR Spectrum Estimation 
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APPENDIX D. COMPUTER SIMULATION CODE 


SESELEETEELSSEEESFESELELEEEEESEEEESEESEEEEEESSESEEFEESES 
%% backfor.m 


%% This program extrapolates the data using forward- 
%% backward linear prediction. The AR spectral 

%% estimation is also calculated. 

&% 
TETEELESEEEEEEESEEEEEEEETEEEEELETLETEELETEEEELETTEETES 


load kmdat.mat; % load data sequence 
load robd.mat; 
load coeff.mat; % load prediction/AR 


% parameters 
=coeff; 
%predf=robd; 
predf=kmdat; 


freql=.1; 
freq2=.2; 
freq3=.21; 


ord=16; % model order 
ext=256; % amount extrapolated 
pad=1024; 


% Forward prediction 


n=65; 

for i=l:ext; 
p(i)=linpred(predf,n,a,ord); 
predf(n)=p(i); 
=n+1; 

end 


% Backward prediction 


flip=flipud(kmdat) ; % extrapolate backwards 
m=65; ° 
for i=l:ext; 
p(i)=linpred(flip,m,a,ord); 
flip(m)=p(i); 
m=m+1 ; 
end 


predb=flipud(flip); 








total=(predb(1:ext) ;predf]; 


clear predf predb flip p a kmdat 


$plot (total) ,title('EXTRAPOLATED SEQUENCE') 
%xlabel('sample'), ylabel('amplitude') 
meta bflp 

tr-use 


k=length (total) 
window=hamming(k) ; 
finham=window .* total; 
clear window total 
pow=fft (finham, pad) ; 
n=length (pow) ; 

pyy=pow .* conj(pow) ./n; 


£=0.5*(O0:n/2)/(n/2); 

pyy (n/2+2:n)=(J; 
pownorm=(1/max(pyy)) .* pyy; 
powlog=10 .* 10g10(pownorm) ; 
clear pow pyy pownorm k 


fl=(freql -100; freqi 100] 
f2=[(freq2 -100; freq2 100] 
f£3=[freq3 -100; freq3 100) 
axis((0 .5 -60 0)); 


~e Me Se 


a=[l;coeff]; 
sar=abs(fft(a,pad)) .°2; 

sar= 1.0 ./ sar}; 
sarnorm=(1/max(sar)) .* sar; 
sarlog= 10 .* 1log10(sarnorm) ; 
sarlog=sarlog(1l:n/2 + 1); 


plot(f,powlog, '-',f1(:,1),f1(:,2),'--',f£2(:,1),£2(:,2)) 


title('PERIODOGRAM OF EXTRAPOLATED SEQUENCE’) 
ylabel('RELATIVE PSD (dB) ') 

Xlabel('FRACTION OF SAMPLING FREQUENCY') 

meta bflp3 

pause 


plot(f,sarlog,'-',f1(:,1),£1(:,2),'-~',f2(:,1),£2(:,2)) 


title('AR SPECTRUM ESTIMATION') 
xlabel('FRACTION OF SAMPLING FREQUENCY ') 
ylabel('RELATIVE PSD (dB) ') 

meta 

pause 

axis; 
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TELSEEFEEEEEESEEEETEEEELEESTEEESEEEELESSESEESGEEETGELSEETEES 
function lp=linpred(x,n,a,ord) 
SEEEEEEEELEEEEEETESTETEELESEESEELTEEEEELEEEEEEEEEEEELEEES 


% linpred computes the nth element of a n-1 sequence using 
% simple linear prediction, x-data sequence; n=element to be 
% predicted; a=filter coefficients; ord=order of predictor. 


i=1; 

for k=1l:ord; 
y(i)=a(k) *x(n-k) ; 
i=itl; 

end 


lp=-sum(y) ; 
SELEEEEEEEEETELELESEEEELELEELETELETETELEEEEEE ELSE LETTS 
%% mfblp.m 

%% 

%% This program computes the prediction/AR parameters 
%% using the modified forward/backward method. 


SETETELEEEEEEETELELEELELELEELETEELEETELTELEELEEL ELE ETE 


S$load kmdat.mat; 
load robd.mat; 


%x=kmdat; 
x=robd; 


n=length (x) ; % number of data points 
m=96; 
np=n-m; 


for i=l:np; % load data matrix A 
for k=1:m; 
A(i,k)=x(itm-k) ; 
A(itnp,k)=x(itk) ; 
end 
end 


%(u,s,Vvj=svda(A); % singular value decomposition 
{eigvec,eigval]=eig(A'*A); 


k=4; % number of sinusoids present 
% in signal, assumed known 


92 





eigval=diag(eigval) ; 
steigval=diag(s) .°2; % eigenvalues of A'A 
w=zeros(m,1); 


for i=l:n-m; 


b(i,1)=x(mti); % 2(n-m)-by-1 desired 
b(itn-m,1)=x(i); % response vector 

end 

for i=1:k; ; 
w=w + (eigvec(:,m-i+1) ./ eigval(m-i+1) *eigvec(:,m-i+1) 

end 

%for i=1:k; 

% =w + (v(:,i) ./ eigval(i) * v(:,1)' * A' * b); 

tend ; 


ferase coeff.mat 

diary coeff.mat 

wW=-W 

diary off 

'q coeff.mat 
TESEETEEELEETEETEETTEELEEETELEEEELEEETEEETEEEELELTEETS 
$% pron.m 

SB 


%% This program finds an IIR filter with a prescribed 
%% time domain response using Prony's Method for filter 
%$% design. The impulse response of that filter is then 
%% computed and a periodogram is computed to estimate 
%% the spectrum of the data. The data is extrapolated 
%% by increasing the impulse duration. 

$% 


SESS T SST TTSTTSSTEELITEETEEEETEETEVELITITEEEEEEEEELEBREE 


load kmdat.mat; 


h=kmdat; % time domain impulse response 
hf=flipud(h); 


nb=16; % order of numerator 
na=16; % order of denominator 
zpad=1024; 

{b,aJ=prony (h,nb,na) ; % IIR filter coefficients 
(d,c]=prony (hf,nb,na) ; 

n=128; 

x=(1 zeros(1,n-1)]; 

y=filter(b,a,x) ; ynew=y (65:n) '; 


yf=filter(d,c,x); ynewf=flipud(yf(65:n)'); 
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